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The present investigation designs a systematic method for finding the latent roots and 
the principal axes of a matrix, without reducing the order of the matrix. It is characterized 
by a wide field of applicability and great accuracy, since the accumulation of rounding errors 
is avoided, through the process of "minimized iterations". Moreover, the method leads to 
a well convergent successive approximation procedure by which the solution of integral 
equations of the Fredholm type and the solution of the eigenvalue problem of linear differ- 
ential and integral operators may be accomplished. 



I. Introduction 

The eigenvalue problem of linear operators is of 
central importance 4 for all vibration problems of 
physics and engineering. The vibrations of elastic 
structures, the flutter problems of aerodynamics, 
the stability problem of electric networks, the 
atomic and molecular vibrations of particle phys- 
ics, are all diverse aspects of the same fundamental 
problem, viz., the principal axis problem of quad- 
ratic forms. 

In view of the central importance of the eigen- 
value problem for so many fields of pint 4 and 
applied mathematics, much thought has been de- 
voted to the designing of efficient methods by 
which the eigenvalues of a given linear operator 
may be found. That linear operator may be of 
the algebraic or of the continuous type; that is, a 
matrix, a differential operator, or a Fredholm 
kernel function. Iteration methods play a prom- 
inent part in these designs, and the literature on 
the iteration of matrices is very extensive. 2 In 
the English literature of recent years the works of 
H. Hotelling [1] 8 and A. C. Aitken [2] deserve 
attention. H. Wayland [3] surveys the field in 
its historical development, up to recent years. 
W. U. Kincaid [4] obtained additional results by 
improving the convergence of some of the classical 
procedures. 

1 The preparation of this paper was sponsored in pari by the office of 
Naval Research. 

2 The basic principles of the various Iteration methods are exhaustively 
treated in the well-known hook on Elementarj matrices by R. A. Frazer, 
\V. J. Duncan, and A. R. Collar (Cambridge University Press, 1938); 
(MacMillan. New York, N. Y., 1947). 

BFigures in brackets Indicate the Literature references at the end of this 
paper. 
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The present investigation, although starting 
out along classical lines, proceeds nevertheless in 
a different direction. The advantages of the 
method here developed 4 can be summarized as 
follows: 

1. The iterations are used in the most economi- 
cal fashion, obtaining an arbitrary number of 
eigenvalues and eigensolutions by one single set 
of iterations, without reducing the order of the 
matrix. 

2. The rapid accumulation of fatal rounding 
errors, common to all iteration processes if applied 
to matrices of high dispersion (large "spread" of 
the eigenvalues), is effectively counteracted by 
the method of "minimized iterations". 

3. The method is directly translatable into 
analytical terms, by replacing summation by 
integration. We then get a rapidly convergent 
analytical iteration process by which the eigen- 
values and eigensolutions of linear differential and 
integral equations may be obtained. 

II. The Two Classical Solutions of 
Fredholm' s Problem 

Since Fredholm's fundamental essay on integral 
equations [5], we can replace the solution of linear 
differential and integral equations by the solution 

4 The literature availahle to the author showed no evidence that the 
methods and results of the present investigation have been found before. 
However, A. M. Ostrowski of the University of Basle and the Institute for 
Numerical Analysis informed the author that his method parallels the , 
earlier work of some Russian scientists; the references given by Ostrowski ' 
are: A. Krylov, Izv. Akad. Nauk SSSR 7, 491 to 539 (1931); N. Luzin, Izv. 
Akad. Nauk SSSR 7, 903 to 958 (1931). On the basis of the reviews of these 
papers in the Zentralblatt, the author believes that the two methods coincide 
only in the point of departure. The author has not, however, read these 
Russian papers. 
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of a set of simultaneous ordinary linear equations 
of infinite order. The problem of Fredholm, if 
formulated in the language of matrices, can be 
stated as follows : Find a solution of the equation 



y—\Ay=b, 



(1) 



where & is a given vector, X a given scalar para- 
meter, and A a given matrix (whose order event- 
ually approaches infinity) ; whereas y is the 
unknown vector. The problem includes the 
inversion of a matrix (X= °°) and the problem of 
the characteristic solutions, also called "eigenso- 
lutions", (b=0) as special cases. 

Two fundamentally different classical solutions 
of this problem are known. The first solution is 
known as the "Liouville-Neumann expansion" 
[6]. We consider A as an algebraic operator and 
obtain formally the following infinite geometric 



series: 



v- 



1 



1-\A 



6=(1+X^+X 2 ^1 2 + . . .)&. Cz) 



This series converges for sufficiently small values 
of |X| but diverges beyond a certain |X| = |Xi|. 
The solution is obtained by a series of successive 
"iterations"; 5 we construct in succession the 
following set of vectors: 

b =b 



b,=Ab 
b 2 =Ab x 



b n +i — Ab n 



(3) 



and then form the sum: 



y = b u + \b i + \% + 



(4) 



The merit of this solution is that it requires 
nothing but a sequence of iterations. The draw- 
back of the solution is that its convergence is 
limited to sufficiently small values of X. 

The second classical solution is known as the 
Schmidt series [7]. We assume that the matrix 
A is "nondefective" (i. e. that all its elementary 
divisors are linear) . We furthermore assume that 

5 Throughout this paper the term "iteration" refers to the application of 
the given matrix A to a given vector b, by forming the product Ab. 



we possess all the eigenvalues 6 n t and eigenvectors 
u t of the matrix A, defined by the equations 



Au i =\x i n i (i=l,2, . . . ,n). 



(5) 



If A is nonsymmetric, we need also the "adjoint" 
eigenvectors u t *, defined with the help of the 
transposed matrix ^4*: 



A*u i *=iiti i * (i=l,2, 
We now form the scalars 
bu t * 



,»). (6) 



Ti B 



U v Uf 



(7) 



and obtain y in form of the following expansion: 



y- 



1 — X/xi 1 — Xju 2 



1— XjUn 



(8) 



This series offers no convergence difficulties, 
since it is a finite expansion in the case of matrices 
of finite order and yields a convergent expansion 
in the case of the infinite matrices associated with 
the kernels of linear differential and integral 
operators. 

The drawback of this solution is — apart from 
the exclusion of defective matrices 7 — that it pre- 
supposes the complete solution of the eigenvalue 
problem associated with the matrix A. 

III. Solution of the Fredholm Problem by 
the S-Expcmsion 

We now develop a new expansion that solves 
the Fredholm problem in similar terms as the 
Liouville-Neumann series but avoids the conver- 
gence difficulty of that solution. 

We first notice that the iterated vectors b , b u 
b 2 , ... cannot be linearly independent of 
each other beyond a certain definite b k . All these 
vectors find their place within the ^-dimensional 
space of the matrix A, hence not more than n of 
them can be linearly independent. We thus know 
in advance that a linear identity of the following 
form must exist between the successive iterations. 



b n +gib m - l +g 2 b m - 2 -\ \-g m b =0. 



(9) 

6 We shall use the term "eigenvalue" for the. numbers m defined by (5), 
whereas the reciprocals of the eigenvalues: \i=llm shall be called "character, 
istic numbers". 

7 The characteristic solutions of defective matrices (i. e. matrices whose 
elementary divisors are not throughout linear) do not include the entire 
^-dimensional space, since such matrices possess less than n independent 
principal axes. 
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We cannot tell in advance what m will be, except 
for the lower and upper bounds : 



1 <*m^n. 



(10) 



How to establish the relation (9) by a systematic 
algorithm will be shown in section VI. For the 
time being we assume that the relation (9) is 
already established. We now define the poly- 
nomial 

G(x)=x m +g l x M - l +- • -+g m , (11) 

together with the "inverted polynomial" (the 
coefficients of which follow the opposite sequence) : 



S«(X)=l+0iX+0 2 X 2 + - • >+g m X" 



(12) 



Furthermore, we introduce the partial sums of the 
latter polynomial: 

S =1 



S 2 (A) = l + <7iX+<7 2 X 2 



S«_i(X) = H-piX+---+ff«-iX" 



y . 



(13) 



We now refer to a formula which can be proved by 
straightforward algebra : s 

+S M _,(X) • xv+. • -+S ■ X"- 1 *— ». 

(14) 

Let us apply this formula operationally, replacing 
x by the matrix A, and operating on the vector 
b . In view of the definition of the vectors &,. 
the relation (9) gives: 



0(A) .6 = 0, 



(15) 



and thus we obtain 

S m (\) 



1— \A 

and hence: 



b Q --=S m7i (\)h+S m ^(\)\b l + --- + S \ m -%„- l , 

(16) 



r 



1 , S m -,(X)& +iS w - 2 (X)X6 1 +- • •+SoX" i -'6 m - 1 



l-\A 



b Q =- 



S m (\) 



(17) 



8 In order to prevent this paper from becoming too lengthy the analytical 
details of the present investigation are kept to a minimum, and in a few places 
the reader is requested to interpolate the missing steps. 



If we compare this solution with the earlier 
solution (4), we notice that the expansion (17) 
may be conceived as a modified form of the Liou- 
ville-Neumann series, because it is composed of 
the same kind of terms, the difference being only 
that we weight the terms X^* by the weight factors 



„„ O w _£_i(\) 

w * — SJST' 



(18) 



instead of taking them all with the uniform weight 
factor 1. This weighting has the beneficial effect 
that the series terminates after m terms, instead of 
going on endlessly. The weight factors w t are 
very near to 1 for small X but become more and 
more important as X increases. The weighting 
makes the series convergent for all values of X. 

The remarkable feature of the expansion (17) is 
its complete generality. No matter how defective 
the matrix A may be, and no matter how the 
vector 6 was chosen, the expansion (17) is always 
valid, provided only that we interpret it properly. 
In particular we have to bear in mind that there 
will always be m polynomials S k (\), even though 
every S k (\) may not be of degree k, due to the 
vanishing or the higher coefficients. For example, 
it could happen that 

G(x)=x m , (19) 

so that 

S m (X) = l + 0X+0X 2 + " . . +0X m , (20) 

£*(X) = 1 + 0X + 0X 2 + . . . +0X", (21) 
and the formula (17) gives: 

y = b + \b l + . . . +\ m ~ i b m . l . (22) 

IV. Solution of the Eigenvalue Problem 

The Liouville-Neumann series cannot give the 
solution of the eigenvalue problem, since the ex- 
pansion becomes divergent as soon as the param- 
eter X reaches the lowest characteristic number 
Xi . The Schmidt series cannot give the solution of 
the eigenvalue problem since it presupposes the 
knowledge of all the eigenvalues and eigenvectors 
of the matrix A. On the other hand, the expan- 
sion (17), which is based purely on iterations and 
yet remains valid for all X, must contain implicitly 
the solution of the principal axis problem. In- 
deed, let us write the right side of (1) in the form 



b=S m (\)b. 



(23) 
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Then the expansion (17) loses its denominator 
and becomes: 

y=S m ^(\)T Q +S m . 2 (\)\b l + . . . + SoX m_1 "^-i. (24) 

We can now answer the question whether a 
solution of the homogeneous equation 



y-\Ay = Q, 



(25) 



is possible without the identical vanishing of y. 
The expression (23) shows that b can vanish only 
under two circumstances; either the vector b, or 
the scalar S m (\) must vanish. Since the former 
possibility leads to an identically vanishing y, 
only the latter possibility is of interest. This 
gives the following condition for the parameter X. 



S m (\)=0. 



(26) 



The roots of this equation give us the character- 
istic values \ = A f , whereas the solution (24) yields 
the characteristic solutions, or eigenvalues, or 
principal axes of the matrix A: 



Ut=S m -i(\i)bo+S m - 2 (\i)\ibi + 



+ SoK 



(27) 



It is a remarkable fact that although the vector 
b was chosen entirely freely, the particular linear 
combination (27) of the iterated vectors has in- 
variant significance, except for an undetermined 
factor of proportionality that remains free, in view 
of the linearity of the defining equation (25). 
That undetermined factor may come out even as 
zero, i. e., a certain axis may not be represented in 
the trial vector b at all. This explains why the 
order of the polynomial S m (\) need not be neces- 
sarily equal to n. The trial vector b may not give 
us all the principal axes of A. What we can say 
with assurance, however, is that all the roots of 
S m (\) are true characteristic values of A, and all 
the Ui obtained by the formula (24) are true char- 
acteristic vectors, even if we did not obtain the 
complete solution of the eigenvalue problem. The 
discrepancy between the order m of the poly- 
nomial G{"^ 
equation 



/x; and the order n of the characteristic 



TO = 



d\\ — M 



&nl 



• &ln 



= 



(28) 



will be the subject of the discussion in the next 
section. 

Instead of substituting into the formula (27), 
we can also obtain the principal axes u t by a 
numerically simpler process, applying synthetic 
division. By synthetic division we generate the 
polynomials : 



X — fli 



-x m - l +g[x m - 2 + . . . +£-1. (29) 



We then replace x 3 ' by bj and obtain: 

u i =b m _ l +g{b m _2+ • • • +gin-ibo. (30) 

The proof follows immediately from the equation 

(A- Mi )u i -G(A)-b = 0. (31) 

V. The Problem of Missing Axes 

Let us assume that we start with an arbitrary 
" trial vector" b and obtain by successive itera- 
tions the sequence: 



bo,b u b 2 , . . . , b n . 



(32) 



Similarly we start with the trial vector b* and 
obtain by iterating with the transposed matrix A* 
the adjoint sequence: 



bo,bl,bl, . . . , K. 



(33) 



Let us now form the following set of "basic 
scalars": 

c i+k =b r bl=b^bt (34) 

The remarkable fact holds that these scalars 
depend only on the sum of the two subscripts i and 
k; e. g. 

b k ^ 1 b k l 1 = bl. 1 b k+1 = b k bt (35) 

This gives a powerful numerical check of the 
iteration scheme since a discrepancy between the 
two sides of (35) (beyond the limits of the round- 
ing errors) would indicate an error in the calcula- 
tion of b k+l or b* k+l if the sequence up to b k and b* k 
has been checked before. 

Let us assume for the sake of the present argu- 
ment that A is a nondefective matrix, and let us 
analyze the vector b in terms of the eigenvectors 
u t , while b* Q will be analyzed in terms of the 
adjoint vectors u*: 



b =T l u 1 J rr 2 u 2 - 
bl = T\u\-\-T 2 u^- 



+ T n U n , 



(36) 
(37) 
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Then the scalars c { become: 



with 



c f =Pijui+P2/4+ • • • +PnfS n , (38) 

Pk = T k T* k . (39) 



The problem of obtaining the p. from the c. is 
the problem of "weighted moments", which can 
be solved as follows: Assuming that none of the 
p k vanish and that all the \ are distinct, we 
establish a linear relation between n-\-l consecu- 
tive c t , of the following form: 



M0 + M1 + . . 

C1V0+C2V1+ • • 



+c„_ 1 7 7n _ 1 + c n =0 
+c n r] n _ 1 + c n+1 = Q 



C n *?0 + Cn+l*h+ . . . +^2»-l^»-l + C 2 » 







y (40) 



Then the definition of the c { shows directly that 
the set (40) demands 



where 



FGO=o, 



(41) 



F(x)=ri + 7i 1 X+ . . . +Vn-iX n ~ 1 + X\ (42) 

Hence, by solving the recurrent set (40) with the 
help of a "progressive algorithm", displayed in 
the next section, we can obtain the coefficients of 
the characteristic polynomial (42), whose roots 
give the eigenvalues p.. 

Under the given restricting conditions none of 
the fit roots has been lost, and we could actually 
establish the full characteristic equation (41). It 
can happen, however, that b is orthogonal to some 
axis u*, and it is equally possible that b* is 
orthogonal to some axis u k . In that case r* and 
r* drop out of the expansions (36) and (37) and 
consequently the expansion (38) lacks both 
pi and p k . This means that the scalars Cj are 
unable to provide all the \x h since p, t and y. k are 
missing. The characteristic equation (41) cannot 
be fully established under these circumstances. 

The deficiency was here caused by an unsuitable 
choice of the vectors b and b* ; it is removable by a 
better choice of the trial vectors. However, we 
can have another situation where the deficiency 
goes deeper and is not removable by any choice of 
the trial vectors. This happens if the fi t roots of 
the characteristic equation are not all distinct. 
The expansion (38) shows that two equal roots 
X f and \ k cannot be separated since they behave 



exactly as one single root with a double amplitude. 
Generally, the weighted moments c t can never 
show whether or not there are multiple roots, 
because the multiple roots behave like single roots. 
Consequently, in the case of multiple eigenvalues 
the linear relation between the c t will not be of the 
n th but of a lower order. If the number of distinct 
roots is m, then the relation (40) will appear in the 
following form: 



Co*?0 + Cl?7l+ • • • +C m -\Vm-l + C m =0 
C1V0 + C2V1+ • • • +Cm7?m-l + Cm+i = 



CmVo + Cm+lVl + 



y .(43) 



+ C 2 m-lVm-l+C 2m = > 



Once more we can establish the polynomial 

G(x)=r 1o + Vl x+ . . . + Vm _ l x m - l +x m , (44) 

but this polynomial is now of only m ih order and 
factors into the in root factors 



(x— m)(x— m) • • • (x—Pm), 



(45) 



where all the /z* are distinct. After obtaining 
all the roots of the polynomial (44) we can now 
construct by synthetic division the polynomials: 



Gix) _ 
x—»k 



1 +glx m - 2 + 



+ £■ 



(46) 



and replacing x 3 by bj we obtain the principal axes 
of both ^1 and A*: 



u k =b m -i+glb m - 2 + 



+gl-A: 



K=K-i+glK-2+ • • • +g k m-ibl 



(47) 



This gives a partial solution of the principal 
axis problem, inasmuch as each multiple root con- 
tributed only one axis. Moreover, we cannot tell 
from our solution which one of the roots is single 
and which one multiple, nor can the degree of 
multiplicity be established. In order to get further 
information, we have to change our trial vectors 
and go through the iteration scheme once more. 
We now substitute into the formulae (47) again 
and can immediately localize all the single roots 
by the fact that the vectors Uj associated with 
these roots do not change (apart from a propor- 
tionality factor), whereas the u k belonging to 
double roots will generally change their direction. 
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A proper linear combination of the new u* k , and the 
previous u k establishes the second axis associated 
with the double eigenvalue n k ; we put 

u\ = u k u 2 k =u k +yu k 

u\*=u* k ul*=u£+y*u* k 

The factors 7 and 7* are determined by the con- 
ditions that the vectors u\ and u 2 k have to be bi- 
orthogonal to the vectors u\* and uf. In the case 
of triple roots a third trial is demanded, and so on. 

An interesting contrast to this behavior of 
multiple roots associated with nondefective matri- 
ces is provided by the behavior of multiple roots 
associated with defective matrices. A defective 
eigenvalue is always a multiple eigenvalue, bub 
here the multiplicity is not caused by the collapse 
of two very near eigenvalues, but by the multi- 
plicity of the elementary divisor. This comes 
into evidence in the polynomial G(x) by giving a 
root factor of higher than first order. Whenever 
the polynomial G(x) reveals a multiple root, we 
can tell in advance that the matrix A is defective 
in these roots, and the multiplicity of the root 
establishes the degree of deficiency. 

It will be revealing to demonstrate these con- 
ditions with the help of a matrix that combines 
all the different types of irregularities that may 
be encountered in working with arbitrary matrices. 
Let us analyze the following matrix of sixth order: 



1 


2 


3 














1 


4 

















1 




















2 













































The eigenvalue 2 is the only regular eigenvalue 
of this matrix. The matrix is "singular", because 
the determinant of the coefficients is zero. This, 
however, is irrelevant from the viewpoint of the 
eigenvalue problem, since the eigenvalue "zero" 
is just as good as any other eigenvalue. More 
important is the fact that the eigenvalue zero is a 
double root of the characteristic equation. The 
remaining three roots of the characteristic equa- 
tion are all 1. This 1 is thus a triple root of the 
characteristic equation; at the same time the 
matrix has a double deficiency in this root, 
because the elementary divisor associated with 



this root is cubic. The matrix possesses only 
four independent principal axes. 

What will the polynomial G(x) become in the 
case of this matrix? The regular eigenvalue 2 
must give the root factor x— 2. The regular 
eigenvalue has the multiplicity 2 but is reduced 
to the single eigenvalue and thus contributes 
the factor x. The deficient eigenvalue 1 has the 
multiplicity 3 but also double defectiveness. 
Hence, it must contribute the root factor (x— l) 3 . 
We can thus predict that the polynomial G(x) 
will come out as follows. 

G(x)=x(x-2)(x—iy=x 5 -5x i +9x*-7x 2 +2x. 

Let us verify this numerically. As a trial vector 
we choose 





b = b 


:=i, 


1,1,1 


,1,1 






The successive 


iterations yield the following: 


b = 


1 


1 


1 


1 


1 


1 


&! = 


6 


5 


1 


2 








b 2 = 


19 


9 


1 


4 








b 3 = 


40 


13 


1 


8 








K= 


69 


17 


1 


16 








b 5 = 


106 


21 


1 


32 








b 6 = 


151 


25 


1 


64 








K= 


1 


1 


1 


1 


1 


1 


K= 


1 


3 


8 


2 








«= 


1 


5 


23 


4 








bl= 


1 


7 


46 


8 








K= 


1 


9 


77 


16 








K= 


1 


11 


116 


32 








K= 


1 


13 


163 


64 









We now construct the c x by dotting b with the 
b* (or 60 with the b t ) ; we continue by dotting b 6 
with b[, . . . bl (or bl with b u . . . b & ) . This gives 
the following set of 2n-\-l = 13 basic scalars. 

Cj = 6, 14, 33, 62, 103, 160, 241, 362, 555, 884, 
1477, 2590, 4735. 
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The application of the progressive algorithm of 
section VI to these b t yields G(x) in the predicted 
form. We now obtain by synthetic divisions: 



Q(x)_„ 4 

x-l 



x i -4x 3 + 5x 2 -2x, 



" r W ^o.4_C 



x — 2 



--x*-3x s + 3x 2 — x, 



^=x*-5x A +9x 2 -7x+2. 



Inverting these 1 polynomials we obtain the matrix 






-2 


5 


-4 


1 








-1 


3 


-3 


1 





2 


-7 


9 


-5 


1 






The product of this matrix with the iteration 
matrix B (omitting the last row b 6 ) yields three 
principal axes Ui] similarly the product of the same 
matrix with the iteration matrix B* yields the 
three adjoint axes u]\ 

u(l) = _8 
u(2) = 2 
u(0)= 2 2 

u*(l)= 0-8 
u*(2)= 2 
u*(0)= 2 2 

Since G(x) is of only fifth order, while the order 
of the characteristic equation is 6, we know that 
one of the axes is still missing. We cannot de- 
cide a priori whether the missing axis is caused by 
the duplicity of the eigenvalue 0, 1, or 2. 9 How- 
ever, a repetition of the iteration with the trial 
vectors 

b =b;=i, i, i, i, i,o, 

causes a change in the row u(0) and u*(Q) only. 
This designates the eigenvalue u=Q &s the double 
root. The process of biorthogonalization finally 

yields: 10 

9 We have in mind the general case and overlook the fact that the over- 
simplified nature of the example makes the decision trivial. 

10 The leader is urged to carry through a similar analysis with the same 
matrix, but changing the 0, diagonal elements of the rows 5 and 6 to l, 0, 
and to 1, 1. 



1^(0) =0 1 1 

tt 2 (0)=0 1 -1 

^(0)=0 1 1 

u 2 (0)=0 1-1 

VI. The Progressive Algorithm for the 
Construction of the Characteristic Poly- 
nomial G(x) 

The crucial point in our discussions was the 
establishment of a linear relation between a 
certain b m and the previous iterated vectors. 
This relation leads to the characteristic poly- 
nominal G(x), whose roots G(m)=Q yield the 
eigenvalues mi- Then by synthetic division we 
can immediately obtain those particular linear 
combinations of the iterated vectors b h which 
give us the eigenvectors (principal axes) of the 
matrix A. 

We do not know in advance in what relation the 
oi-dcr m of the polynomial G(x) will be to the order 
n of the matrix A. Accidental deficiencies of the 
trial vectors b 07 b* , and the presence of multiple 
eigenvalues in A can diminish m to any value 
between 1 and ii . For this reason we will follow 
a systematic procedure that generates G(x) 
gradually, going through all degrees from 1 to m. 
The procedure comes automatically to a halt 
when the proper m has been reached. 

Our final goal is to solve the recurrent set of 
equations: 



C i7o + CiT?i + 



+ c m =0 



C m riQ-{-C m+ i7ii+ . . . +C 2 



(48) 



This is only possible if the determinant of this 
homogeneous set vanishes: 



C Ci . 
C\ c 2 . 



t'TO + l 



Cm 
Cm+l 



. C 2n 



=o. 



(49) 
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Before reaching this goal, however, we can cer- 
tainly solve for any k<^m the following inhomog- 
eneous set (the upper k is meant as a superscript) : 






C k Vo + C k + 1 Vi- 



+ c k = 
+ c k+l = 



V • 



+ c 2k =h k > 



(50) 



The freedom of h k in the last equation removes the 
overdetermination of the set (48). The proper 
m will be reached as soon as h m turns out to be 
zero. 

Now a recurrent set of equations has certain 
algebraic properties that are not shared by other 
linear systems. In particular, there exists a 
recurrence relation between the solutions of three 
consecutive sets of the type (50). This greatly 
facilitates the method of solution, through the 
application of a systematic recurrence scheme 
that will now be developed. 

We consider the system (50) and assume that 
we possess the solution up to a definite k. Then 
we wi]l show how this solution may be utilized 
for the construction of the next solution, which 
belongs to the order k-\-\. 

Our scheme becomes greatly simplified if we 
consider an additional set of equations that 
omits the first equation of (50) but adds one 
more equation at the end: 



Citi+c 2 rj k 2 -\ r-c* + i = 



(51) 



Let us now multiply the set (50) by the factor 



l^k + l 



(52) 



and add the set (51). We get a new set of equa- 
tions that can be written down in the form: 



CoV k o + 1 + C l vi +1 + - 



•+c k+1 = Q 



C k Vo + 1 + C^iV k i + 1 +'-+C 2 ^ 



(53) 



provided that we put 



^ +1 = 



= <lkVo 



--q k v k i+vi 



(54) 



ni +1 =&-i+«J 

We now evaluate the scalar 

Ck+\r]o +1 -\ \-C 2k+ ir) k k + 1 + C 2k + 2 = h k + l (55) 

which is added to (53) as the last equation. 

WTiat we have accomplished is that a proper 
linear combination of the solutions 17* and rfi 
provided us with the next solution r?f +1 . But 
now exactly the same procedure can be utilized 
to obtain rj$ +1 on the basis of rf* and qj +1 . 

For this purpose we multiply the set (51) by 

k 



a*=- 



flk+l 



(56) 



and add the set (53), completed by (55) but 
omitting the first equation. 
This gives: 



cA w +Cie +1 +---+c* 







^+it?i +1 +c,-k^ +1 H \-c 2k+2 =0. 

provided that we put: 



(57) 



5?r=2,^+^ +l ^ 



t2 + 1 =l k t2 + vV 



-k+i — pT 
Vk+i — <lk 



(58) 



Once more we evaluate the scalar 



c t+ ,tf +1 +- 



'~T^2k + 3 % 



(59) 



which is added to (57) as the last equation. 

This analytical procedure can be translated into 
an elegant geometrical arrangement that generates 
the successive solutions iyf and 77? in successive 
columns. The resulting algorithm is best ex- 
plained with the help of a numerical example. 

For this purpose we choose the eigenvalue prob- 
lem of an intentionally over-simplified matrix, 
since our aim is not to show the power of the 
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method but the nature of the algorithm, which 
leads to the establishment of the characteristic 
equation. The limitation of the method due to 
the accumulation of rounding errors will be dis- 
cussed in the next section. 
Let the given matrix be: 




We iterate with the trial vector b =l, 0, 0, and 
obtain: 

1 

13 4 7 

28 24 12 

208 64 112 



We transpose the matrix and iterate with the 
trial vector b * = l, 0, 0, obtaining: 

1 

13 5 -23 

28 -4 -20 

208 80 -368 

We dot the first row and the last row with the 
opposing matrix and obtain the basic scalars c% as 
follows: u 

1, 13, 28, 208, 448, 3328, 7168. 

These numbers are written down in a column, 
and the scheme displayed below is obtained. 








0.5 


1 


1.5 


2 


2.5 


3 


hi = 


1 


13 


-141 


147. 6923075 


163. KM2569 








Q' = 

1 


-13 
1 


10. 84615384 


1.047-103173 


1.106382990 









13 




1 


-13 










28 






1 


-2. 15384616 


-13.617021249 






208 








1 


-1. L06382987 


-16.000000004 





448 










1 


0.000000003 


-16 


3328 












1 





7 1 f 58 














1 



(60) 



Instead of distinguishing between the 77, and rji 
solutions we use a uniform procedure but mark 
the successive columns alternately as "£1111" and 
"half columns"; thus we number the successive 
columns as zero, one-half, one, . . . The scheme 
has to end at Sifull column, and the end is marked 
by the vanishing of the corresponding "head- 
number" h t . In our scheme the head-number is 
zero already at the half-column 2.5, but the scheme 
cannot end here, and thus we continue to the 
column 3, whose head-number becomes once more 
0, and then the scheme is finished. The last 
column gives the polynomial G(x), starting with 
the diagonal term and proceeding upward: 

Gtx) = l.x* + 0.x 2 -16x+0, 
=x 3 — 16x. 

The head-numbers //, are always obtained by 
dotting the column below with the basic column 
d) e. g., at the head of the column 2 we find the 



number -163.4042569. This number was ob- 
tained by the following cumulative multiplication: 

448-1 +208- (-1.106382987)4-28. (-13. 617021 249). 

The numbers '/, represent the negative ratio of 
two consecutive A, numbers: 



2r=" 



"hV 



11 Instead of iterating with A and A* n times, we can also iterate with A 
alone In times. Any of the columns of the iteration matrix can now be 
chosen as ci numbers since these columns correspond to a dotting of the 
iteration matrix with bl = l, 0, 0, . . ., respectively 0, 1, 0, 0, . . .; 0, 0, 1, 0, 

0, . . .; and so on. The transposed matrix is not used here at all. E. C. 
Bower of the Douglas Aircraft Co. points out to the author that from the 
machine viewpoint a uniform iteration scheme of 2n iterations is preferable 
to a divided scheme of n+n iterations. The divided scheme has the advan- 
tage of less accumulation of rounding errors and more powerful checks on the 
successive iterations. The uniform scheme hus the advantage that more 
than one column is at our disposal. Accidental deficiencies of the b* Q vector 
can thus be eliminated, by repeating the algorithm with a different column. 
(For this purpose it is of advantage to start with the trial vector 6o=l, 1, 

1, . . . 1.) In the case of a symmetric matrix it is evident that after n itera- 
tions the basic scalars should be formed, instead of continuing with n more 
iterations. 



904428 — 50- 
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e - g., 2i.5=1.106382990 was obtained by the fol- 
lowing division: 

-(-163.4042569) 
147.6923075 

The scheme grows as follows. As soon as a 
certain column C% is completed, we evaluate the 
associated head-number h t ; this provides us with 
the previous ^-number g 2 _ Ka . We now construct 
the next column C i+H by the following operation. 
We multiply the column (7 f _ H by the constant 
gi-x an d add the column d: 

However, the result of this operation is shifted 
down by one element; e. g. in constructing the 
column 2.5 the result of the operation 



1.10638299- (-2.15384616) + 

(-13.617021249) = - 



-16, 



is not put in the row where the operation occurred, 
but shifted down to the next row below. 

The unfilled spaces of the scheme are all 



zero 



The outstanding feature of this algorithm is that 
it can never come to premature griej, provided only 
that the first two c-numbers, c and c h are different 
from zero. Division by zero cannot occur since 
the scheme comes to an end anyway as soon as the 
head-number zero appears in one of the full 
columns. 

Also of interest is the fact that the products of 
the head-numbers associated with the full columns 
give us the successive recurrent determinants of 
the c t ) e.g., the determinants 



1 



and 



1 


13 




1 


13 


28 


13 


28 




13 


28 


208 








.28 


208 


448 


1 


13 


28 


208 




13 


28 


208 


448 




28 


208 


448 


3328 




208 


448 


33 


28 


7168 





are given by the successive products 

1, l-(— 141) = — 141, (-141H- 163.4042569) = 
23040, and 23040-0=0, 



13 


28 




13 


28 


208 


28 


208 




28 


208 


448 








208 


448 


3328 



Similarly the products of the head-numbers of 
the half-columns give us similar determinants, but 
omitting c from the sequence of onumbers. In 
the example above the determinants 



13, 



are given by the products: 

13, 13-147.6923075 = 1920, 1920-0 = 

The purpose of the algorithm (60) was to gener- 
ate the coefficients of the basic identity that exists 
between the iterated vectors b t . This identit} r 
finds expression in the vanishing of the poly- 
nomial G(x): 



G(x)=0 



(61) 



The roots of this algebraic equation give us the 
eigenvalues of the matrix A. In our example we 
get the cubic equation 

x 3 — 16x=0, 

which has the three roots 

Ml = 0, M2 =4, M3=— 4. (62) 

These are the eigenvalues of our matrix. In order 
to obtain the associated eigenvectors, we divide 
G(x) by the root factors: 



G(x) 



--x 2 -W 



<?(*). 



=x 2 +4x 



#(*). 



-4 " ' *~ X+4: 

This gives, replacing x k by b k : 

^(0) = -166 +&2, 

u(-4) = ~46 1 +6 2 . 
Consequently, if the matrix 



-4:X. 



( 



-16 





1 





4 


1 





— 4 


1 



is multiplied by the matrix of the b t (omitting 6 3 ), 
we obtain the three eigenvectors Ui\ 
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16 





1\ 


/ 1 





°\ 


/ 12 


24 





4 


1 


13 


4 


7 = 


( 80 


40 





-4 


1/ 


\28 


24 


12/ 


\-24 


8 



12\ /tt(0) \ 
40 = «(4) 
-16/ V(-4)/. 



If the same matrix is multiplied by the matrix of the transposed iterations bt* (omitting 63), we obtain 
the three adjoint eigenvectors u\: 




l) (13 5 -23 W 

1/ \28 -4 -20/ \- 



12 


-4 


- 20\ 


/u*(p) 


80 


16 


-112 1 = 


( i/*(4) 


24 


-24 


72/ 


V*(-4) 



The solution of the entire eigenvalue problem is thus accomplished. 



VII. The Method of Minimized Iterations 

In principle, the previous discussions give a com- 
plete solution of the eigenvalue problem. We have 
found a systematic algorithm for the generation of 
the characteristic polynomial O(fi). The roots of 
this polynomial gave the eigenvalues of the matrix 
A. Then the process of synthetic division estab- 
lished the associated eigenvectors. Accidental de- 
ficiencies were possible but could be eliminated by 
additional trials. 

As a matter of fact, however, the "progressive 
algorithm" of the last section has its serious limi- 
tations if large matrices are involved. Let us 
assume that there is considerable "dispersion" 
among the eigenvalues, which means that the ratio 
of the largest to the smallest eigenvalue is fairly 
large. Then the successive iterations will grossly 
increase the gap, and after a few it era I ions the 
small eigenvalues will be practically drowned out. 
Let us assume, e. g., that we have a 12 by 12 mat- 
rix, which requires 12 iterations for the generation 
of the characteristic equation. The relatively mild 
ratio of 10:1 as the "spread" of the eigenvalues is 
after 12 iterations increased to the ratio 10 12 :1, 
which means that we can never get through with 
the iteration scheme because the rounding errors 
make all iterations beyond the eighth entirely 
valueless. 

As an actual example, taken from a physical 
situation, let us consider four eigenvalues, which 
are distributed as follows: 



1, 



50, 2000. 



Let us assume, furthermore, that we start with a 
trial vector that contains the four eigenvectors in 
the ratio of the eigenvalues, i. e., the eigenvalue 
2000 dominates with the amplitude 2000, compared 
with the amplitude of the eigenvalue 1. After 
one iteration the amplitude ratio is increased to 
4-10 6 , after two iterations to 8-10 9 . The later 



iterations can give us no new information, since 
they practically repeat the second iteration, multi- 
plied every time by the factor 2000. The small 
eigenvalues 1 and 5 are practically obliterated and 
cannot be rescued, except by an excessive accuracy 
that is far beyond the limitations of the customary 
digital machines. 

We will now develop a modification of the cus- 
tomary iteration technique that dispenses with 
this difficulty. The modified scheme eliminates 
the rapid accumulation of rounding errors, which 
under ordinary circumstances destroys the value 
of high order iterations. The new technique pre- 
vents the large eigenvalues from monopolizing the 
scene. It protects the small eigenvalues by con- 
stantly balancing the distribution of amplitudes in 
the most equitable fashion. 

As an illustrative example, let us apply this 
method of "minimized iterations" to the above- 
mentioned dispersion problem. If the largest 
amplitude is normalized to 1, then the initial 
distribution of amplitudes is characterized as 
follows: 

0.0005, 0.0025, 0.025, 1. 

Now, while an ordinary iteration would make 
this distribution still more extreme, the method 
of minimized iterations changes the distribution 
of amplitudes as follows: 

0.0205 0.1023 1 -0.0253. 

We see that it is now the third eigenvector that 
gets a large weight factor, whereas the fourth 
eigenvector is almost completely in the back- 
ground. 

A repetition of the scheme brings about the 
following new distribution: 

0.2184 1 -0.1068 0.0000. 
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It is now the second eigenvector that gets the 
strongest emphasis. 

The next repetition yields: 



I 



-0.2181 0.0018 0.0000, 



and we see that the weight shifted over to the 
smallest eigenvalue. 

After giving a chance to each eigenvalue, the 
scheme is exhausted, since we have all the infor- 
mation we need. Consequently the next mini- 
mized iteration yields an identical vanishing of 
the next vector, thus bringing the scheme to its 
natural conclusion. 

In order to expose the principle of minimized 
iterations, let us first consider the case of symmetric 
matrices: 

A*=A. (63) 

Moreover, let us agree that the multiplication of 
a vector b by the matrix A shall be denoted by a 
prime: 

Ab = V. (64) 

Now our aim is to establish a linear identity 
between the iterated vectors. We cannot expect 
that this identity will come into being right from 
the beginning. Yet we can approach this identity 
right from the beginning by choosing that linear 
combination of the iterated vectors b' Q and b 
which makes the amplitude of the new vector as 
small as possible. Hence, we want to choose as 
our new vector &i the following combination: 



bi = K~ a b 0l 
where a Q is determined by the condition that 



This gives 



(b' — a &o) 2 = minimum. 



b' o b 

"0 



Notice that 



bi-b =0, 



(65) 

(66) 
(67) 

(68) 



i.e. the new vector b x is orthogonal to the original 
vector b . 

We now continue our process. From bi we 
proceed to b 2 by choosing the linear combination 

b 2 = b[-aA-^A, (69) 

and once more a x and /3 are determined by the 



condition that b\ shall become as small as possible. 
This gives 






(70) 



A good check of the iteration b[ is provided by the 
condition 



b , X=b l b' =b\. 



(71) 



Hence, the numerator of p Q has to agree with the 
denominator of a!. 

The new vector b 2 is orthogonal to both b and bi. 

This scheme can obviously be continued. The 
most remarkable feature of this successive mini- 
mization process is, however, that the best linear 
combination never includes more than three terms. 
If we form 6 3 , we would think that we should put 



h = b2—a 2 b 2 — pibi — yobo 



(72) 



But actually, in view of the orthogonality of b 2 to 
the previous vectors, we get 



_b' 2 b _b 2 b' 
To — T2~— ~W~ u - 



K 



(73) 



Hence, every new step of the minimization process 
requires only two correction terms. 

By this process a succession of orthogonal vectors 
is generated: 12 



bo, &i, b 2 , . . . , b„ 



(74) 



unti] the identity relation becomes exact, which 
means that 



b m =0. 



(75) 



If the matrix A is not symmetric, then we 
modify our procedure as follows: We operate 
simultaneously with A and A*. The operations 
are the same as before, with the only difference 
that the dot products are always formed between 
two opposing vectors. The scheme is indicated 
as follows: 



b 
bi = b f — a o b Q 



bl = b* ' — a b* 



^bobp^ bobo 
a ° b b* b* b 



12 The idea of the successive orthogonalization of a set of vectors was prob- 
ably first employed by O. Szasz, in connection with a determinant theorem 
of Hadamard; cf. Math, es phys. lapok (in Hungarian) 19, 221 to 227, (1910). 
The method found later numerous important applications. 
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b 2 = b , 1 — a l b 1 — p {) b 



b^bl' — onbl—Pobo 



_b[Vj_b*ibi 

ai ~b 1 bi~b:b l 

a b[b* bl'bo b^bi 



(76) 



b Q bl b* b b* o b 
b2 = b2—a 2 b 2 —fi]bi bl = bl. — a 2 bl—fiib* f etc. 

Operationally the prime indicates multiplica- 
tion by the matrix A. Hence, the succession of 
bi vectors represents in fact a successive set of 
polynomials. Replacing A by the more familiar 
letter x, we have: 

& =3 'b Q 

bi=(x—ao)b 



b 2 =(x—ai)b l —Pob 

bz=(x — a 2 )b 2 — Pibi 



(77) 



b m = (x — a m -i)b m _i — p m _ 2 b m _ 2 =0 



This gradual generation of the characteristic 
polynomial G{x) is in complete harmony with the 
procedure of the "progressive algorithm", dis- 
cussed in section VI. In fact, the successive poly- 
nomials of the set (77) are identical with the 
polynomials found in the full columns of the progres- 
sive algorithm (60). This explains the existence of 
the recurrence relation 

Pm+l(x) = (x — a n )Pn(x)— ft»-l2>«-l(a0, (78) 

without additional 7, 5, . . . terms. The 
existence of such a relation is a characteristic 
feature of the recurrent set of equations that are at 
the basis of the entire development. 

Although the new scheme goes basically through 
the same steps as the previously discussed "pro- 
gressive algorithm", it is in an incomparably 
stronger position concerning rounding errors. 
Apart from the fact that the rounding errors do 
not accumulate, we can effectively counteract their 
influence by constantly checking the mutual 
orthogonality of the gradually evolving vectors 
bi and b t . Any lack of orthogonality, caused by 
rounding errors, can immediately be corrected by 
the addition of a small correction term. 13 By this 
procedure the orthogonality of the generated 



vector system does not come gradually out of gear, 
However, quite apart from the numerical advan- 
tages, the biorthogonality of the vectors b t and b] 
has further appeal because it imitates the behavior 
of the principal axes. This is analytically an 
eminently valuable fact that makes the transition 
from the iterated vectors to the principal axes a 
simple and strongly convergent process. 

In order to see the method in actual operation, 
let us apply it to the simple example of section 
VI. Here the matrix A is of third order, and 
thus we have to construct the vectors b ,b h b 2 , 63, 
and the corresponding adjoint vectors. We ob- 
tain the following results: 





6 = 1 





K= 1 







K= 13 4 


7 


«' = 13 


5 -23 


&i= 


4 


7 


K= 


5 -23 


b[ = 


-141 -28 


-79 


«'= 


141 -69 279 




1677 
* 1= -T41 




ft 


--; 41 - hi 



--11.89361702 
b 2 = 0, 19.57446808, 4.25531914 
&;=0, -9.53191490, 5.44680854 
&; = 0, -17.02127656, 3.40425542 
b* 2 '=0, 16.34042562, -32.68085142 



«2 



180.78768715 
-163.4042553 

-1.106382981 

63 = 0, 0, 



ft= 



-163.40425746 



-141 
= 1.158895443 
«=0, 0, 







13 See section IX, eq 90. 



The associated polynomials become: 

Pi(x)=x— 13 

p 2 (aO = (x+ 11.89361702) (x--13) + 141 

=x 2 - 1.10638298^- 13. 61702126 
p z (x) = (x + 1.106382981) (x 2 + 10638298*- 
13.61702126) — 1.158895443 (x— 13) 

=* 3 -16x 
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Comparison with the scheme (60) shows that 
the coefficients of these very same polynomials 
appear in the full columns 0, 1, 2, 3, of the pro- 
gressive algorithm of section VI. 

VIII. Solution of the Eigenvalue Problem 
by the Method of Minimized Iterations 

The biorthogonal property of the vector system 
b t , b* leads to an explicit solution of the eigen- 
value problem, in terms of the vectors b t . Let us 
first assume that the matrix A is of the noneffec- 
tive type, and let us analyze the vectors b t in terms 
of the eigenvectors u t . The method by which the 
vectors b t were generated, yields directly the 
relation 

b i =Pi(ni)u l +p i (ti 2 )u 2 + . . . +p t (nm)y> m . (79) 

If this relation is dotted with u* k , we obtain, in 
view of the mutual orthogonality of the two sets 
of axes: 

bi-u* k =pi(n k )u k -ul (80) 

Let us now reverse the process and expand the 
u i in terms of the b t : 



/ Ui = a i ob -{-aiibi-\- . . . -{- a im -\b m _i. 

The dotting by b* k yields: 

u r bl 



ctik~- 



(81) 



(82) 



00 <?l <?2 



+ Pm-l(Hi)- 



bk'K 

Let us denote the "norm" of b k by <r k : 

<rjc=b k >b* k , (83) 

while the norm of u k will be left arbitrary. Then 
the expansion (81) becomes: 

(84) 

This expansion contains the solution of the prin- 
cipal axis problem. The eigenvectors Ui are gen- 
erated in terms of the vectors b t , which are the 
successive vectors of the process of minimized 
iterations. The expansion (84) takes the place of 
the previous "^-expansion" (27), which solved the 
eigenvector problem in terms of the customary 
process of iteration. 

The adjoint axes are obtained in identical 
fashion: 



(Tq a i c 2 



-p m -l(Hi) 



(85) 



The expansion (84) remains valid even in the 
case of defective matrices. The only difference is 
that the number of principal axes becomes less 
than n since a multiple root p j} if substituted into 
(84) and (85) cannot contribute more than one 
principal axis. 14 However, a defective matrix ac- 
tually possesses less than n pairs of principal axes, 
and the above expansions give the general solution 
of the problem. 

An interesting alternative of the expansion (84) 
arises if we go back to the original Fredholm 
problem and request a solution in terms of the 
minimized vectors b if rather than the simply iter- 
ated vectors of the expansion (17). One method 
would be to make use of the Schmidt series (8), 
expressing the u t of that series in terms of the b u 
according to the expansion (84). However, the 
Schmidt series holds for nondefective matrices 
only, while we know that a solution must exist 
for any kind of matrices. 

Hence, we prefer to proceed in a somewhat 
different fashion. We expand y directly in terms 
of the vectors b t : 



2/=2/ &o+yi6i + . • .+y m -ib m . 1 . 



(86) 



We substitute this expansion into the equation 
(1), replacing b' k by 



&*=&*+i + aA + ftfc-i&* 



(87) 



Then we compare coefficients on both sides of the 
equation. The result can be described as follows: 
Let us reverse the sequence of the a r coefficients 
and let us do the same with the j8 -coefficients. 
Hence, we define 

a =a m -i 

(88) 



~-a 



i* The reader is urged to carry through the process of minimized iterations 
and evaluation of the principal axes for the defective matrix 



(: i :) 

\0 0/ 



which has only one pair of principal axes; (choose the trial vector in the form 
6o=&; = l, 1, 1). 
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We now construct the following "reversed" set 
of polynomials: 

Pi(x) =X—a {) 

p 2 (x) = (x-aOpiiz) -ft, (89) 



= G(x) 

Then the solution of the Frcdholm problem (1) 
is given by the following expansion : 

V^~QT\ [bm-i+Pi(v)b m -2(u)b m - 2 + • • • +Pm-i(u)bo], 



where we have put 



(90) 

(91) 



The expansion (90) is completely general and 
remains valid, no matter how the vector b of the 
right side was given, and how regular or irregular 
the matrix A may be. The only condition to be 
satisfied is that the vector b* — while otherwise 
chosen arbitrarily — shall be free of accidental 
deficiencies, i. e., b* shall not be orthogonal to 
some u k if b is not simultaneously orthogonal to 

The expansion (90) leads once more to a solu- 
tion of the eigenvector problem, this time obtained 
with the help of the "reversed" polynomials 

p t (x): 



u i =b m - l J rPi{iXx)bm-2+ 



+P»-i(m<)&o. (92) 



The expansions (92) and (84) actually coincide— 
except for a factor of proportionality — due to 
algebraic reasons. 




1.10638298 

5.10638298 

-2.89361702 




In order to see a numerical example for this 
solution of the eigenvalue problem let us return 
once more to the simple problem of section VI. 
The minimized b t and b* vectors associated with 
this matrix were given at the end of section VII. 
together with the associated polynomials p t (x). 
We now construct the reversed polynomials 
Pi(x). For this purpose we tabulate the a t 
and (3 t : 

13 -141 

-11.89361702 1.158895443 
-1.106382981 

We reverse the sequence of this tabulation: 



-1.106382981 
-11.89361702 
13 



1.58895443 
-141 



and construct in succession: 

Po=l 

p x (x) =x+ 1.106382981 

/J 2 (a0 = (x+1 1.8361702)/J, (*) — 1.58895443 

=.r 2 +13;r+12 
p z (x) = (x--lZ)p 2 (x) + Ul 

=x 3 — 16x 

The last polynomial is identical with p 3 (x)—G(x). 
The zeros of this polynomial arc: 

Mi = 0, M2 = 4, /x 3 =— 4; 

substituting these values into ^(m), P\{v) > Po we 
obtain the matrix: 



12 1.10638298 

80 5.10638298 

, — 24 -2.89361702 



■) 



The product of this matrix with the matrix of the 
b f vectors gives the three principal axes u t : 





4 

19.57446808 



\ / 12 24 

7 1=1 80 40 

4.25531914/ \-24 8 



12\ M0) 

40 1=1 u (4) 
-16/ V(-4), 
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in complete agreement with the previous result 
of section VI, but now obtained by an entirely 
different method. If the 6-matrix is replaced by 
the &*-matrix, the adjoint axes ^*(0), u*(4). 
u*( — 4) are obtained. 

IX. The Lateral Vibrations of a Ear 

In order to study the power of the method in 
connection with a vibration problem of large 
dispersion, the elastic vibrations of a bar were 
investigated. The bar was clamped at one end 
and free at the other. Moreover, the cross section 
of the bar changes in the middle, as shown in 
figure 1. 




Figure 1. Elastic bar, clamped at one end, free at the 
other; the moment of inertia makes a jump in the middle 
in the ratio 1:2. 

The change of the cross section was such that the 
moment of inertia jumped from the value 1 to 2. 
The differential equation that describes the vibra- 
tions of such a bar is the following fourth-order 
equation: 



d 2 
dx" 



[m d2 y~\ 



dx 2 J 



-\y, 



(93) 



with the boundary conditions: 

y(0)=0 y"(l)=o 
y'(0)=0 y"'(l)=0 



and 



k(x) = \(o^x<~\k(x)=2(^<Cx^l\ 



(94) 



(95) 



The differential operator d/dx was replaced by the 
difference operator A/ Ax, with Ai = l. The length 
of the bar was chosen as 1=13, thus leading to a 
12 by 12 matrix; (since y(0)=y(l) = 0). 

The first step was the inversion of the matrix. 
This was easily accomplished, since a matrix that 
is composed of a narrow band around the diagonal 
can be inverted with little labor. The eigenvalues 
Hi cf the inverted matrix are the reciprocals of the 
original a*: 

Py = ~ (96) 

A; 



The general theory has shown that the iteration 
scheme applied to an arbitrary matrix automati- 
cally yields a biorthogonal set of vectors b t and b*; 
they can be conceived as the building blocks from 
which the entire set of principal axes may be gener- 
ated. In the present problem dissipative forces 
are absent, which makes the matrix A symmetric 
and the problem self-adjoint. Hence, 



bi = b*, 



(97) 



and we get through with a single set of iterations. 

Now the general procedure would demand that 
we go through 12 minimized iterations before the 
stage &i2=0 is attained. However, the study of a 
system with high dispersion has shown that in 
such a system the method of minimized iterations 
practically separates the various vibrational 
modes, starting with the highest eigenvalue and 
descending systematically to the lower eigenvalues, 
provided that we employ a trial vector b which 
weights the eigenvectors approximately accord- 
ing to the associated eigenvalues. In the present 
problem the trial vector 1, 0, . . ., was not used 
directly but iterated with the matrix A, and then 
iterated again. The thus obtained vector b" 
was employed as the b of the minimized iteration 
scheme. 

The strong grading of the successive eigen- 
vectors has the consequence that in k minimized 
iterations essentially only the highest k vibrational 
modes will come into evidence. This is of eminent 
practical value, since it allows us to dispense with 
the calculation of the very low eigenvalues (i. e., 
very high frequencies since we speak of the eigen- 
values of the inverted matrix), which are often of 
little physical interest, and also of little mathe- 
matical interest, in view of the fact that the 
replacing of the d operator by the A operator be- 
comes in the realm of high frequencies more and 
more damaging. 

Whether the isolation actually takes place or 
not can be tested with the help of the p t (x) poly- 
nomials, which accompany the iteration scheme. 
The order of these polynomials constantly increases 
by 1. The correct eigenvalues of the matrix A 
are obtained by evaluating the zeros of the last 
polynomial p m (x)=0. What actually happens, 
however, is that the zeros of the polynomials 
Pi(x) do not change much from the beginning. 
If the dispersion is strong, then each new poly- 
nomial basically adds one more root but corrects 
the higher roots by only small amounts. It is 
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thus well possible that the series of largest roots 
in which we are primarily interested is practically 
established with sufficient accuracy after a few 
iterations. Then we can stop, since the later 
iterations will change the obtained values by negli- 
gible amounts. The same can be said about the 
vibrational modes associated with these roots. 

This consideration suggests the following succes- 
sive procedure for the approximate determination 
of the eigenvalues and eigenvectors (vibrational 
modes) of a matrix. As the minimization scheme 
proceeds and we constantly obtain newer and 
newer polynomials pi(x), we handle the last poly- 
nominal obtained as if it were the final polynomial 
p m (x). We evaluate the roots of this polynomial 
and compare them with the previous roots. 
Those roots that change by negligible amounts are 
already in their final form. 

A similar procedure holds for the evaluation of 
the eigenvectors u t . Here the biorthogonality of 
the vectors b t and 6* — which is reduced to simple 
orthogonality in the case of a symmetric matrix — 
is of very great help. Let us assume that I lie 
length of the vectors b t is normalized to 1, by 
replacing b t by bj-yjeji. Then the expansions (84) 
and (85) show that the following matrix must be 
an orthogonal — although in the diagonal terms not 
normalized — matrix : 



1 




£2.0*1) 

V0-2 


Pm-lM 
V°"m-1 


1 


V°"i 


^2(M2) 

Vo-2 


Pm-l 0*2) 

V°m-i 



s.v<ro 



£i0*») 



p2(Pm) 



I <*2 



Pm-lOfm) 



(98) 



The dot-product of any two rows of this matrix 
must come out as zero — thus providing us with a 
powerful check on the construction of the p t poly- 
nomials and the correctness of the roots fi i} which 
are the roots of the equation p m (n)=Q. 15 In the 
case of strong dispersion, the transformation 
matrix (98) is essentially reduced to the diagonal 
terms and one term to the right and to the left of 
the diagonal; i. e., the eigenvector u k is essentially 
a linear combination of three 6-vectors only, viz., 
h- 2 , 6*-i. and b k . 



These general conditions are well demonstrated 
by the tabulation of the final results of the above- 
mentioned bar problem. The minimized itera- 
tions were carried through up to m~6. On the 
basis of these iterations the first six eigenvalues 
and the first five vibrational modes of the clamped- 
free bar were evaluated. The iterations were con- 
stantly watched for orthogonality. After obtain- 
ing a certain b iy this b t was immediately dotted 
with all the previous bj. If a certain dot-product 
b t - bj came out as noticeably different from zero, 
the correction term 



bj-bj h 
" bj* h 



(99) 



was added to b u thus compensating for the in- 
fluence of rounding errors. By this procedure the 
10 significant figure accuracy of the calculations 
was constantly maintained. 18 

The roots of the successive polynomials p t (x) are 
tabulated as follows: 



Mi 


M2 


M3 


2256. 926071 
. 943939 
. 943939 
. 943939 
. 943939 
. 943939 






48. 1610705 
. 20 57755 
. 2037825 
. 2037825 
. 2037825 


"T272§Il428 

. 355958260 
. 356269794 

. 3562699S0 


M4 


M5 


Me 




















1.513923859 

. 582259337 
. 5829955952 






0. 546327303 
.591117817 




0. 2498132719 



15 Duo to reasons of algebra, the orthogonality of the matrix (98) holds not 
only for the final m hut for any value of m. 

16 In a control experiment that imitated the conditions of the vihrating bar 
but with a more regular matrix, the results were analytically predictable and 
the computational results open to an exact check. This example vividly 
demonstrated the astonishing degree of noninterference of orthogonal vectors. 
The spread of the eigenvalues was 1 : 3200. The trial vector 6o strongly 
over-emphasized the largest eigenvalue, containing t he lowest and t he highest 
eigenvectors with an amplitude ratio of 1 : 10 8 ; (this means that if the vector 
of the smallest eigenvalue were drawn with the length of 1 in., the vector of 
the largest eigenvalue, perpendicular to the previous one, would span the 
distance from Los Angeles to Chicago). The slightest inclination between the 
two vectors would fatally injure the chances of the smallest eigenvalue. 
When the entire computation scheme was finished, the analytically required 
eigenvalues were computed and the comparison made. The entire set of 
eigenvalues, including the last, agreed with a maximum error of 2 units in the 
niiilh significant figure, thus demonstrating that the method is completely 
free of the cumulation of rounding errors. The author is indebted to 
Miss Fannie M. Gordon, of the National Bureau of Standards, for the 
eminently careful and skilled performance of the computing operations. 
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The successive orthogonal transformation ma- 
trices (98) show likewise strong convergence. 
We tabulate here only the last computed trans- 
formation matrix (rounded off to four decimal 
places), which expresses the first six eigenvectors 

r 1 
-.0028 





We notice how quickly the elements fall off to 
zero as soon as we are beyond one element to the 
right and one to the left of the main diagonal. 
The orthogonal reference system of the b t and 
the orthogonal reference system of the u t are 
thus in close vicinity of each other. 

The obtained five vibrational modes U\ . . ., u 5 
(u& being omitted since the lack of the neigh- 
bour on the right side of the diagonal makes the 
approximation unreliable) are plotted graphically 
in figure 2. 



0028 











1 


.0316 


.0004 





0316 


1 


. 1497 


.0081 


0044 


-. 1520 


1 


.2693 


0010 


.0335 


-. 2793 


1 


0002 


-.0087 


.0779 


-.3816 




2 4 6 8 10 12 14 

Figure 2. The first five lateral vibrational modes of the elastic 
bar, shown in figure 1 . 

mi=2256.94; ju 2 =48.20; ^3=5.36; ju4=1.58; ms=0.59 

X. The Eigenvalue Problem of Linear In- 
tegral Operators 

The methods and results of the past sections 
can now be applied to the realm of continuous 
operators. The kernel of an integral equation 
can be conceived as a matrix of infinite order, 
which may be approximated to any degree of 
accuracy by a matrix of high but finite order. 
One method of treating an integral equation is 
to replace it by an ordinary matrix equation of 



Ui, . . ., Ui in terms of the first six nor- 
malized bij-yjo-i vectors, making use of the roots 
of p 6 (u)=0. The diagonal elements are nor- 
malized to 1 : 



' 




.0249 
.4033 

1 



sufficiently high order. This procedure is, from 
the numerical standpoint, frequently the most 
satisfactory one. However, we can design meth- 
ods for the solution of integral equations that 
obtain the solution by purely analytical tools, 
on the basis of an infinite convergent expansion, 
such as the Schmidt series, for example. The 
method we are going to discuss belongs to the 
latter type. We will find an expansion that is 
based on the same kind of iterative integrations 
as the Liouville-Neumann series, but avoiding the 
convergence difficulties of that expansion. The 
expansion we are going to develop converges 
under all circumstances and gives the solution of 
any Fredholm type of integral equation, no 
matter how defective the kernel of that integral 
equation may be. 17 

Let us first go back to our earlier method of 
solving the Fredholm problem. The solution 
was obtained in form of the ^-expansion (17). 
The difficulty with this solution is that it is based 
on the linear identity that can be established 
between the iterated vectors b t . That identity 
is generally of the order n\ if n grows to infinity, 
we have to obtain an identity of infinite order 
before our solution can be constructed. That, 
however, cannot be done without the proper 
adjustments. 

The later attempt, based on the method of 
minimized iterations, employs more adequate 
principles. We have seen that to any matrix, 
A, a biorthogonal set of vectors b t and b\ can 
be constructed by successive minimizations. The 
set is uniquely determined as soon as the first 
trial vectors b and b* are given. In the case 
of the inhomogeneous equation (1) the right side 
b may be chosen as the trial vector b , while b 
is still arbitrary. 

i? The Volterra type of integral equations, which have no eigenvalues and 
eigensolutions, are thus included in our general considerations. 
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The construction of these two sets of vectors 
is quite independent of the order n of the matrix. 
If the matrix becomes an integral operator, the 
bi and 6* vectors are transformed into a bior- 
thogonal set of functions : 



</> (x), </>i (x), <h(x), . . 

</>o», <t>l(x), <&(x), . . 



(100) 



which are generally present in infinite number. 
The process of minimized iterations assigns to 
any integral operator such a set, after cj> (x) and 
<f> Q (x) have been chosen. 

Another important feature of the process of 
minimized iterations was the appearance of a 
successive set of polynomials #<(/*), tied together 
by the recurrence relation 

P <+1 (m) = (m— ai)p<(/i)— Pi-iPi-i(ii). (101) 

This is again entirely independent of the order n 
of the matrix A and remains true even if the 
matrix A is replaced by a Frcdholm kernel K(x^). 
We can now proceed as follows. We stop at 
an arbitrary p m (x) and form the reversed set of 
polynomials p«Gi), defined by the process (88). 
Then we construct the expansion: 



ym(x) = 






[0m-l(aO+2>l(M)0ro-2(z)+ • 



(102) 



This gives a successive approximation process that 
converges well to the solution of the Fredholm 
integral equation 



y(x)—\Ky(x)=<l>o(x). 



In other words: 



y(x) = \imy m (x). 



(103) 



(104) 



By the same token we can obtain all the eigen- 
values and eigcnsolutions of the kernel K(x£), 
if such solutions exist. For this purpose we 
obtain the roots /z* of the polynomial p m (n) by 
solving the algebraic equation 



p m (n)=0. 



(105) 



The exact eigenvalues /x f of the integral operator 
K(x,g) are obtained by the limit process: 



\imp m (nt)=0, 



(106) 



where the largest root is called fx l} and the sub- 
sequent roots are arranged according to their 
absolute values. The corresponding eigenfunc- 
tions are given by the following infinite expansion: 



(x) = lim [- 



.(*) 



o-o 



-Pifai) 



nt(x) = ] 

n—Kn |__ °" 

, . </) m -i(x) "| 
Cfm-l J 



i(x) 



<*\ 



(107) 



where Mi is the i th root of the polynomial p m (ij). ls 
As a trial function <t> (x) we may choose for 
example 



</>o = const. = l. 



(108) 



However, the convergence is greatly speeded up 
if we first apply the operator K to this function, 
and possibly iterate even once more. In other 
words, we should choose fo=K*l, or even (j) —K 2 'l 
as the basic trial function of the expansion (107). 

We consider two particularly interesting ex- 
amples that are well able to illustrate the nature 
of the successive approximation process here dis- 
cussed. 

(A) The vibrating membrane. In the problem 
of the vibrating membrane we encounter the fol- 
lowing self-adjoint differential operator: 

-fcWh (109) 

This leads to the eigenvalue problem 

-^(av')=X0(O3*£l), (HO) 

with the boundary condition 

y(i)=o. (in) 

The solution of the differential equation (110) is 



y=J d (2^}<x), 



(112) 



where J (x) is the Bessel-f unction of the order zero. 
The boundary condition (111) requires that X shall 
be chosen as follows: 



x,=f 



where £ f are the zeros of J Q (x). 



(IU) 



18 This expansion is not of the nature of the Neumann series, because the 
coefficients of the expansion are not rigid but constantly changing with the 
number m of approximating terms. 
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Now the Green's function of the differential 
equation (110) changes the differential operator 
(109) to the inverse operator, which is an integral 
operator of the nature of a symmetric Fredholm 
kernel function K(x, £). Our problem is to obtain 
the eigenvalues and eigenf unctions of this kernel. 

If we start out with the function <l> (x) = l, the 
operation K<t> gives 

1— x, 

and repeating the operation we obtain 



4' 



^ 3 



and so on. The successive iterations will be 
polynomials in x. Now the minimized iterations 
are merely some linear combinations of the ordi- 
nary iterations. Hence the orthogonal sequence 
<j>i{x) will become a sequence of polynomials of 
constantly increasing order, starting with the 
constant <£o=l- This singles out the 4> k (x) as the 
Legendre polynomials P k (x), but normalized to 
the range to 1, instead of the customary range 
— 1 to +1. The ^normalization of the range 
transforms the polynomials P k (x) into Jacobi 
polynomials G k (p, q; x), with p = g=l [8], which 
again are special cases of the Gaussian hypergeo- 
metric series F(a, /3, y; x), in the sense of F(k+\, 
—k, 1 ; x) ; hence, we get: 



4>o 



= 1 



<j )l (x) = l—2x 

<t> 2 (x) = l-6x + 6x 2 

03 ( x ) = 1 _ 1 2x + 3 Ox 2 - 20r* 



(114) 



The associated polynomials ptix) can be obtained 
on the basis of the relation: 

K<l> m = <l>n+l-\~<Xn<t>n-\~fin-\<l>n-l- (H5) 

This gives: 

l>i{x)=2x — 1 

(116) 
p 2 (x)=24x 2 - 18x+l 

p. 6 (x)=720x 3 -600x 2 +72x-\ J 



In order to obtain the general recurrence relation 
for these polynomials it is preferable to follow the 
example of the algorithm (60) and introduce the 
"half-columns" in addition to the full columns. 
Hence, we define a second set of polynomials q\ k (x) 
and set up the following recurrence relations: 



p n (x)=nxq n _i(x)— p n -i(x) 
Zn(x)=2(2n+l)p n (x)—q n - 1 (x) 



(117) 



We thus obtain, starting with p =l and %=2, 
and using successive recursions: 



Po=l 

Pi(x)=2x— 1 

p 2 (x) = 2±x 2 -l8x+l 

p^(x)=720x 3 - 

600rM-72x-l 



q =2 

q l (x) = 12x—S 

q 2 (x) = 24:0x 2 — 192x + 18 

q*(x) = 10080x 3 -8640x 2 + 
1200^-32 



. . . (118) 

The zeros of the p m (x) polynomials converge to 
the eigenvalues of our problem, but the conver- 
gence is somewhat slow since the original function 
0o=l does not satisfy the boundary conditions 
and thus does not suppress sufficiently the eigen- 
functions of high order. The zeros of the q t poly- 
nomials give quicker convergence. They are 
tabulated as follows, going up to q 5 (x) : 



(119) 



The last row contains the correct values of the 
eigenvalues, computed on the basis of (113): 



Ml 


M2 


M3 


M4 


M5 


0. 6667 
.69155 
. 69166016 
. 6916602716 
. 6916602760 










0. 1084 
. 130242 
. 1312564 
.13127115 








. 035241 
.051130 
. 0532911 






. 014842 

. 025582 


"." 00729 


. 6916602761 


. 13127123 


. 0534138 


. 028769 


.01794 



Mi r 



"£* 2 



(120) 



We notice the excellent convergence of the scheme. 

The question of the eigenfunctions of our prob- 
lem will not be discussed here. 

(B) The vibrating string. Even modes. Another 
interesting example is provided by the vibrating 
string. The differential operator here is 
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dx 2 ' 
with the boundary conditions 
?/(±D=0. 
The solution of the differential equation 

under the given boundary conditions is: 



(121) 



(122) 



(123) 



y f =cos (2i-\-l) -^ x (even modes) . (124) 

y 3 =sin jwx (odd modes). (125) 

This gives the eigenvalues: 



X ,=(?i±i,rY (even modes), (12 



and 



A,= (jx) 2 



<;> 



(odd modes). (127) 



If we start with the trial function $0=1, we will 
get all the even vibrational modes of the string, 
while <I>q=x will give all the odd vibrational modes. 
We first choose the previous alternative. 

Successive iterations give 



r 2 1 
2 2 

24 4^24. 



(128) 



and we notice that the minimized iterations will 
now become a sequence of even polynomials. The 
transformation x 2 = £ shows that these polynomials 
are again Jacobi polynomials G k (p, q; x 2 ), but now 
P = Q = h and we obtain the hyp ergeome trie func- 
tions F{k J r \,— k, \)x 2 ): 



*o=l 

4> l (x) = l-- 3x 2 
<fe(2)=3— 30:z 2 +35;z 4 

3 (x)=5-lO5x 2 +315.z 4 -231x 6 



y (129) 



Once more we can establish the associated 
polynomials p t (x), and the recurrence relation by 
which they can be generated. In the present case 
the recurrence relations come out as follows : 



(130) 



p n (x) = (4:n— l)xq n -i(x)— p n -\{x) 
q_n(x) = (4n+l)p n (x)—q n ^(x) 
starting with p = 1, q = 1 . This yields 

^o=l 
Pi(x)=Sx— 1 
p 2 (x) = \05x 2 -4:5x+l 
p 3 (x) = 10395x 3 -4725x 2 -}-210x-l 

q l (x) = l5x — Q 

q 2 (x)=945x 2 -420x+l[y 

g 3 (z) = 135135r 3 -62370z 2 +3150.r-28 



The successive zeros of the q t (x) polynomials, 
up to q 5 (x), arc tabulated as follows: 
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Mi 


M 2 


M 3 


M 4 


M 5 


0. 40000 
. 405059 
. 405284733 
. 4052847346 
. 4052847346 










6. 02351 

UIIS.V, 

. 04503010 
.0450311)322 








0.011397 
. 015722 
. 016192 






0. 00455 
. 00752 


~6.~662i6~~ 


. 4052847346 


. 0450316371 


.016211 


. 00827 


1 III.-,! HI 



The last row contains the correct eigenvalues, 
calculated from the formula: 



u, 



V>+1 7r/ 



(132) 



The convergence is again very conspicuous. 

(C) The vibrating string. Odd modes. In the 
case of the odd modes of the vibrating string, the 
orthogonal functions of the minimized itera- 
tions are again related to the Jacobi polynomials 
Gk(P,QL;x), but now with _p = <Z=3/2. Expressed in 
terms of the hypergeometric series we now get the 
polynomials of odd orders 



xF\k + 2> -k>2 : x )' 
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0O = £ 

<j) 1 (x)=3x — 5x 3 

(I> 2 (x) = l5x-70x 3 + 63x 5 

3 (x)=35x-315x 3 +693x 5 -429x 7 



► (133) 



The associated Pi(x) polynomials are generated 
by the following recurrence relations: 



(134) 



p n (x) = (4n+ l)xg„_i(aO —p„-.ifr) 
q n (x) = (4:n+S)p n (x) — q n -i(x) 
starting with p = 1, #0=3. We thus get 

Po=l 

Pi(x) = 15x — 1 

p 2 (x)=94:5x 2 -105x'\-l 

^ 3 fe)-135135^ 3 -17325x 2 +378x-l 

go=3 
q l (x) = 105x— 10 

g 2 (x) = 10395a: 2 - 1260.T + 21 
g 3 (x)=2027025x 3 -270270x 2 +6930x-36 



The table of the zeros of g<(#), up to q 5 (x), is 
given as follows: 
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Mi 


M 2 


M 3 


M 4 


M 5 


0. 0952 
.10126 
. 10132106 
.1013211836 
.1013211836 










0. 01995 
. 02500 
. 025323 
. 02533024 








0. 00701 
. 01068 
.011215 






0. 00307 
. 00550 




. 00156 


. 1013211836 


. 0253302P6 


.011258 


. 00633 


. 00405 



The last row contains the correct eigenvalues 
calculated on the basis of the formula: 



1 



(136) 



XL The Eigenvalue Problem of Linear Dif- 
ferential Operators 

Let Dy(x) be a given linear differential operator, 
with given homogeneous boundary conditions of 



sufficient number to establish an eigenvalue prob- 
lem. The problem of finding the eigenvalues and 
eigenf unctions of this operator is equivalent to the 
problem of the previous section in which the eigen- 
value problem of linear integral operators was in- 
vestigated. Let us assume that we know the 
Green's function K(x,£) of the differential equation 



Dy= P . 



(137) 



Then K is the reciprocal operator of D which pos- 
sesses the same eigenfunctions (principal axes) as 
the operator D, whereas the eigenvalues of K are 
the reciprocals of the eigenvalues of D. 

Hence, in principle the eigenvalue problem of 
differential operators needs no special investiga- 
tion. Actually, however, the situation in most 
cases is far less simple. The assumption that we 
are in the possession of the Green's function asso- 
ciated with the differential equation (137) is often 
of purely theoretical significance. Even very 
simple differential operators have Green's func- 
tions, which are outside the limits of our analytical 
possibilities. Moreover, even if we do possess 
the integral operator K in closed form, it is still 
possible that the successive integrations needed 
for the construction of the successive orthogonal 
functions <t>\(x), </> 2 W, ... go beyond our analyt- 
ical facilities. 

In view of this situation we ask the question 
whether we could not relax some of the practically 
too stringent demands of the general theory. We 
may lose somewhat in accuracy, but we may gam 
tremendously in analytical operations if we can 
replace some of the demands of the general theory 
by more simplified demands. The present section 
will show how that may actually be accomplished. 

Leaving aside the method of minimized itera- 
tions, which was merely an additional tool in our 
general program, the basic principle of our entire 
investigation, if shaped to the realm of integral 
operators, may be formulated as follows: 

We start out with a function f (x), which may 
be chosen as/ (aO = l. We then form by iterated 
integrations a set of new functions 

f 1 (x)=Kf (x)J 2 (x)=Kf 1 (x), . . . J m (x)=Kf m _ l (x). 

(138) 

Then we try to establish an approximate linear 
relation between these functions, as accurately as 
possible. For this purpose we made use of the 
method of least squares. 

We notice that the general principle involves 
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two processes: (a) the construction of the iterated 
set (138); (b) the establishment of a, close linear 
relation between them. It is the first process 
where the knowledge of the integral operator 
K=D~ 1 is demanded. But let us observe that 
the relation between the successive /^-functions 
can be stated in reverse order. We then get: 



f m (X) J m _i (x) = Df m (X) , 



,Mx)=Dfi(x). (139) 



If we start with the function f m (x) f then the suc- 
cessive functions of lower order can be formed 
with the help of the given D operator and we can 
completely dispense with the use of the Green's 
function. 

Now the freedom of choosing f (x) makes also 
f m {x) to some extent <i free function. Yet, the 
successive functions f t (x) do not have the same 
degree of freedom. Although f (x) need not sat- 
isfy the given boundary conditions, fi(x) of neces- 
sity satisfies these conditions, whereas/ 2 (#) satisfies 
them more strongly, since not only / 2 (x) but even 
Df 2 (x) satisfies the given boundary conditions. 
Generally, we can say that an arbitrary J„(x) need 
not satisfy any definite differential or integral 
equation, but it is very restricted in the matter of 
boundary conditions: it has to satisfy the bound- 
ary conditions "to the n in order." This means 
that not only j n (x) itself, but the whole series of 
functions: 

j n (x)<Df n (x),D%(x), . . . D^Ux), (140) 

have to satisfy the given boundary conditions. 

To construct a function f n (x) with this property 
is not too difficult. We expand j n (x) into a linear 
set of powers, or periodic functions, or any other 
kind of functions we may find adequate to the 
given problem. The coefficients of this expansion 
will be determined by the boundary condi- 
tions that are satisfied by f n (x) and the iterated 
functions (140). This leads to the solution of 
linear equations. In fact, this process can be 
systematized to a regular recurrence scheme that 
avoids the accumulation of simultaneous linear 
equations, replacing them by a set of separated 
equations, each one involving but one unknown. 

We have thus constructed our set (138), al- 
though in reverse order. We did not use any 
integrations, only the repeated application of the 
given differential operator D. The first phase of 
our problem is accomplished. 



We now turn to the second phase of our program, 
viz., the establishment of an approximate linear 
relation between the iterated functions fi(x). The 
method of least squares is once more at our dis- 
posal. However, here again we might encounter 
the difficulty that the definite integrals demanded 
for the evaluation of the a t and 0* are practically 
beyond our means. Once more we can simplify 
our task. The situation is similar to that of 
evaluating the coefficients of a Fourier series. 
The "best" coefficients, obtained by the method 
of least squares, demand the evaluation of a set 
of definite integrals. Yet we can get a practically 
equally close approximation of a function /(x) by 
a finite trigonometric series f(x), if we use the 
method of "trigonometric interpolation." In- 
stead of minimizing the mean square of the error 
j(x)—j{x), we make/(x) equal toj(x) at a sufficient 
number of equidistant points. This leads not to 
integrations but to simple summations. 

The present situation is quite analogous. To 
establish a linear relation between the/f(a?) means 
that the last function/*^) is to be approximated 
by a linear combination of the previous functions. 
Instead of using the method of least squares for 
this approximation we can use the much simpler 
method of interpolation, by establishing a linear 
relation between the successive j t (x) in as many 
equidistant points as we have coefficients at our 
disposal. For the sake of better convergence, it 
is preferable to omit/ (#) — which does not satisfy 
i he boundary conditions and thus contains the 
high vibrational modes too pronouncedly — and 
establish the linear relation only from /i(x) on. 
For example, if we constructed a trial function 
jz(x) which, together with the iterated/ 2 (x) =Df^(x) 
and doubly iterated fi(x)=D 2 f(x) satisfies the 
given boundary conditions, then we can choose 
two points of the region, e.g., the two endpoints, 
where a linear relation of the form 



Ux)+aUx)+^ l (x)=0, 



(141) 



shall hold. This gives the characteristic poly- 
nomial 6(x) in the form 



G(x)=X 2 +aX + p. 



(142) 



The two roots of this polynomial give us an ap- 
proximate evaluation of the two highest n u (or 
the two lowest X<), i.e., /ii=lAi> an d M£— I/X2. 
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whereas the corresponding eigensolutions are 
obtained by synthetic division: 



G(x) 



X — pi 

G(x) 



X — H2 



==g[x+g f 2 
g'i'x+9'2 



(143) 



which gives 



u l (x)=g[f 2 (x)+g'J } (x) 
u 2 {x) =g"f 2 (x) +g f 2 f fi(x) 



(144) 



(The last root and its eigenfunction are always 
considerably in error, and gives only rough in- 
dications.) 

The remarkable feature of this method is that 
it completely avoids any integrations, requiring only 
the solution of a relatively small number of linear 
equations. 

The following application of the method demon- 
strates its practical usefulness. The method was 
applied to obtain the first three eigenvalues of the 
lateral vibrations of a uniform bar, clamped at 
both ends. The given differential operator is here : 



Dy 



_d?y 

~dx* 



with the boundary conditions 

y(±l)=0 y'(±l)=0. 

Only the even modes were considered, expanding 
y(x) into even powers of x. The approximations 
were carried out to first, second, and third order. 
The obtained eigenvalues are tabulated in the 
following table, the last row containing the correct 
eigenvalues given by Rayleigh [9] 



Ml 


M2 


M3 


0. 0323413 
. 0319686 
. 031963958 






0. 0007932 
. 0010875 




0. 0000788 


0. 031963996 


0. 0010946 


0. 0001795 



We notice that the general convergence behavior 
of this method is exactly the same as that of the 
analytically more advanced, but practically much 
more cumbersome, method of minimized iterations. 



XII. Differential Equations of Second 
Order; Milne's Method 

If a linear differential equation of second order 
with two-end boundary conditions is changed into 
a difference equation and then handled as a matrix 
problem, singularly favorable conditions exist for 
the solution of the eigenvalue problem. The ma- 
trix of the corresponding difference equation con- 
tains only diagonal terms plus one term to the right 
and one to the left. If we now start to iterate with 
the trial vector 



6o=l,0,0, 



0, 



(145) 



we observe that the successive iterations grow by 
one element only, as indicated in the following 
scheme where the dots stand for the nonvanishing 
components : 



bo 

bo = 



(146) 



b n -l = 

K = 



Under these conditions the establishment of the 
linear identity between the iterated vectors is 
greatly simplified since it is available by a succes- 
sive recurrence scheme. The coefficients of the 
equation 



b n + gib n -l+g2bn-2 + 



+g n b =0, (147) 



are directly at our disposal, since the last column 
of the last two vectors gives g lt then the previous 
column gives g 2 , . . . , until finally the first 
column gives g n . The construction of the basic 
polynomial G(x) is thus accomplished and the 
eigenvalues X z 19 directly available by finding the 
roots of the equation 6r(X)=0. 

Professor W. E. Milne of the Oregon State Col- 
lege and the Institute for Numerical Analysis, 
applied the general theory to this problem, but 
with the following modification. Instead of iterat- 

19 We call the eigenvalues here X,-, since the operator D is not inverted into 
an integral operator K. 
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ing with the given matrix A, Milne considers the 
regular vibration problem 



& 2 u , ^ n 



(148) 



where the operator I) has the following signifi- 



cance : A 



Du= [£ 2 + p(x )^ + g {x )~] 



u(x,t). (149) 



The differential equation (148) is now converted 
into a difference equation, with Ax=At = h. 
Then the value of u(ih y jh) are determined by 
successive recursions, starting from the initial 
conditions 



and 



u(ih, 0) = 1, 0, 0, 0, . . . , 0, 



(150) 



u(ih, h)=u(ih f -h). (151) 

The linear identity between the n-\-l vectors 

u(ih, 0), u(ih, A), . , . ,u(i.h,nh), (152) 

leads to a trigonometric equation for the character- 
istic frequencies v if of the following form 



cos nyJi^An^i cos (n—l)vih-\- 
Wc then put 

\i = vl 



+^o=0. 

(153) 

(154) 



On the other hand, the regular iteration method 

gives the eigenvalues X, of the operator Au, de- 
fined in harmony with the operator Du but with 
the modification that the operation d/dx is re- 
placed by the operation A/ Ax. The X- are in the 
following relation to the v i of eq 153 and 154: 



\,=I?= 



sin £ Vi 



(155) 



It is of interest to see that the values (154) of 
Milne are much closer to the true eigenvalues than 
the values obtained by iterations. The values 
of Milne remain good even for high frequencies, 
whereas the iteration method gives gradually 
worse results; this is to be expected since the error 
committed by changing the differential equation 

2<J See J. Research NBS 45, 245 1950 RP2132. The variable "j" of Milne is 
changed to x and his X 2 to X, to avoid conflicts with the notations of the 
present paper. 



to a difference equation must come into evidence 
with ever increasing force, as we proceed to the 
vibrational modes of higher order. 

The following table illustrates the situation. 
It contains the results of one of Milne's examples; 
("Example 1"). Here 

'£_ 

,dx 2 



— G 



•b> 



(156) 



with the boundary conditions: 
Moreover, h was chosen as „ 



(157) 
and n = 7. The 



column y\ k gives the correct frequencies, the 
column V^I gives the frequencies obtained by 
Milne's method, while the column -\\ k gives the 
frequencies obtained by the iteration method. 



k 

1 


V x * 


VXfc 


v** 


:;. 2969 


3. 2898 


3. 2667 


2 


(i. 3623 


6. 3457 


•; |so6 


3 


9. 4777 


9. 45(17 


8.9107 


4 


12.6061 


12. 5664 


11.3138 


5 


15.7398 


15.6820 


13 2S'.I| 


6 


18.8761 


18. 7870 


14.7581 


7 


22. 0139 


21.8430 


15. 6629 



Actually, it is purely a matter of computational 
preference whether we follow the one or the other 
sell erne since there is a rigid relation between the 
two schemes. The frequencies v i obtained by 
Milne's method are in the following relation to 
the frequencies v t obtained by the matrix iteration 
method : 



. h 
sin 2 Vi 



2 Vi 



-v t 



(158) 



Hence, the results obtained by the one scheme can 
be translated into the results of the other scheme, 
and vice versa. 

This raises the question, why it is so beneficial 
to transform the frequencies v t of the Au operator 
to the frequencies v i by the condition 



sin ^ v t -- 



li- 



(159) 



The answer is contained in the fact that the cor- 
rection factor 



. h 
sin 7> Vi 



2 Vi 



(160) 
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is exactly the factor that compensates for the 
transition from du/dx to An/ Ax, if u(x) is of the 
following form : 



u (x) = Ci sin (y t x + Oi) , 



(161) 



where the constants d and 6 t are arbitrary. 

Now it so happens that for large frequencies v i 
the first term of the operator (158) strongly over- 
shadows the other terms. The differential equa- 
tion of the eigenvalue problem for large v i thus 
becomes asymptotically: 



dx 2 



+ v t 2 u t =0, 



(162) 



the solution of which is given by (161). This 
asymptotic behavior of the solution for large 
frequencies makes it possible to counteract the 
damaging influence of the error caused by the 
initial transition to the difference equation. The 
correction is implicitly included in Milne's solu- 
tion, while the results of the matrix iteration 
scheme can be corrected by solving eq 159 for 
the v, 21 

XIII. More-Dimensional Problems 

The present investigation was devoted to 
differential and integral operators that belonged 
to a definite finite range of the variable x. This 
variable covered a one-dimensional manifold of 
points. However, in many problems of physics 
and engineering the domain of the independent 
variable is more than one-dimensional. A few 
general remarks may be in order as to the 
possibility of extending the principles and methods 
of the present investigation to manifolds of higher 
dimensions. 

Although the general theory of integral equa- 
tions reveals that the fundamental properties of 
an integral equation are essentially independent 
of the dimensionality of the variable x, yet from 
the practical viewpoint the eigenvalue problem 
of higher dimensional manifolds does lead to 
difficulties that are not encountered in manifolds 
of one single dimension. The basic difference is 
that an essentially more-dimensional manifold of 



2i This experience is valuable, since in many eigenvalue problems similar 
conditions hold; the eigenfunctions of large order can often be asymptotically 
estimated, in which case the error of the A-process may be effectively cor- 
rected. For example the values m found in section IX for the lateral vibra- 
tions of an inhomogeneous bar may be corrected as follows: 

Hi uncorrected: 2256.944, 48.2038, 5.3563, 1.5830, 0.59, (0.25) 

Hi corrected: 2258.924, 48.4977, 5.4577, 1.6407, 0.62, (0.28) 



eigenvalues is projected on a one-dimensional 
manifold, thus causing a strong overlapping of 
basically different vibrational modes. A good 
example is provided by the vibrational modes 
of a rectangular membrane. The eigenvalues are 
here given by the equation 

\=ai 2 mi 2 + a2 2 m 2 2 , 

where m 1 and m 2 are two independent integers? 
while a x and a 2 are two constants determined by 
the length and width of the membrane. 

As another illustration, consider the bewildering 
variety of spectral terms that can be found within 
a very narrow band of frequencies, if the vibra- 
tional modes of an atom or a molecule are studied. 
To separate all these vibrational modes from each 
other poses a difficult problem that has no analogue 
in systems of 1 degree of freedom where the differ- 
ent vibrational states usually belong to well sep- 
arated frequencies. 

It is practically impossible to expect that one 
single trial function shall be sufficient for the 
separation of all these vibrational states. Nor 
does such an expectation correspond to the actual 
physical situation. The tremendous variety of 
atomic states is not excited by one single exciting 
function but by a rapid succession of an infinite 
variety of exciting functions, distributed accord- 
ing to some statistical probability laws. To imi- 
tate this situation mathematically means that we 
have to operate with a great variety of trial 
functions before we can hope to entangle the very 
dense family of vibrational states associated with 
a more than one-dimensional manifold. 

In this connection it seems appropriate to say 
a word about the physical significance of the 
"trial function' ' cj> (x) that we have employed for 
the generation of an entire system of eigenfunc- 
tions. At first sight this trial function may appear 
as a purely mathematical quantity that has no 
analogue in the physical world. The homogeneous 
integral equation that defines the eigenvalues, and 
the eigenfunctions, of a given integral operator, 
does not come physically into evidence since in 
the domain of physical reality there is always a 
"driving force" that provides the right side of the 
integral equation; it is thus the inhomogeneous 
and not the homogeneous equation that has direct 
physical significance. 

If we carefully analyze the method of successive 
approximations by which the eigenvalues and the 
eigenfunctions of a given integral operator were 
obtained, we cannot fail to observe that we have 
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basically operated with the inhomogeneous equa- 
tion (103) and our trial function ^(x) serves 
merely as the "exciting function" or ' 'driving force." 
Indeed, the solution (107) for the eigenfunctions 
is nothing but a special case of the general solu- 
tion (102), but applied to such values of the param- 
eter X, which make the denominator zero. This 
means that we artificially generate the state of 
"resonance", which singles out one definite eigen- 
value \t and its associated eigenf unction Qi^x). 

From this point of view we can say that, while 
the separation of all the eigenfunctions of a multi- 
dimensional operator might be a practically 
insuperable task — except if the technique of 
"separation" is applicable, which reduces the more- 
dimensional problem to a succession of one- 
dimensional problems — yet it might not be too 
difficult to obtain the sold lion of a given more- 
dimensional integral ('(inn! ion if (he right side 
(i. e., physically, the "driving force") is given as 
a sufficiently smooth function that does not con- 
tain a too-large variety of eigenfunctions. Then 
the convergence of the method may still suffice 
for a solution that gives the output function with 
a practically satisfactory accuracy. This is the 
situation in many antenna and wave-guide 
problems that are actually input-output problems, 
rather than strict resonance problems. In other 
words, what we want to get is a certain mixture 
of weighted eigenfunctions, which appear physi- 
cally together, on account of the exciting mecha- 
nism, while the isolation of each eigenfunction 
for itself is not demanded. Problems of this 
type are much more amenable to a solution than 
problems that demand a strict separation of the 
infinite variety of eigenfunctions associated with 
a multi-dimensional differential or integral oper- 
ator. To show the applicability of the method 
to problems of this nature will be the task of a 
future investigation. 

XIV. Summary 

The present investigation establishes a syste- 
matic procedure for the evaluation of the latent 
roots and principal axes of a matrix, without 
constant reductions of the order of the matrix. 
A systematic algorithm (called " progressive 
algorithm") is developed, which obtains the linear 
identity between the iterated vectors in succes- 
sive steps by means of recursions. The accuracy 



of the obtained relation increases constantly, 
until in the end full accuracy is obtained. 

This procedure is then modified to the method 
of "minimized iterations", in order to avoid the 
accumulation of rounding errors. Great accuracy 
is thus obtainable even in the case of matrices 
that exhibit a large dispersion of the eigenvalues. 
Moreover, the good convergence of the method in 
the case of large dispersion makes it possible to 
operate with a small number of iterations, obtain- 
ing m successive eigenvalues and principal axes 
by only m+1 iterations. 

These results are independent of the order of 
the matrix and can thus be immediately applied 
to the realm of differential and integral operators. 
This results in a well-convergent approximation 
method by which the solution of an integral equa- 
tion of the Fredholm type is obtained by succes- 
sive iterations. The same procedure obtains the 
eigenvalues and eigensolutions of the given integral 
operator, if these eigensolutions exist. 

In the case of differential operators, the too- 
stringent demands of the least square method may 
be relaxed. The approximate linear identity be- 
tween the iterated functions may be established 
by interpolation, thus dispensing with the evalua- 
tion of definite integrals. Moreover, the itera- 
tions may be carried out with the given differen- 
tial operator itself, instead of reverting to the 
Green's function, which is frequently not available 
in closed form. The entire procedure is then free 
of integrations and requires only the solution of 
linear equations. 



The present investigation contains the results 
of years of research in the fields of network analy- 
sis, flutter problems, vibration of antennas, solu- 
tion of systems of linear equations, encountered 
by the author in his consulting and research work 
for the Boeing Airplane Co., Seattle, Wash. The 
final conclusions were reached since the author's 
stay with the Institute for Numerical Analysis, 
of the National Bureau of Standards. The 
author expresses his heartfelt thanks to C. K. 
Stedman, head of the Physical Research Unit of 
the Boeing Airplane Co. and to J. H. Curtiss, 
Acting Director of the Institute for Numerical 
Analysis, for the generous support of his scientific 
endeavors. 
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